/**                   GNU GENERAL PUBLIC LICENSE
 *                       Version 3, 29 June 2007
 *
 * Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
 * Everyone is permitted to copy and distribute verbatim copies
 * of this license document, but changing it is not allowed.
 * See License.txt file that comes with this distribution
 *
 *
 * Author: Vitaly German <vitaly.german@gmail.com>, (C) 2006-2008
 * 
 **/
#include <algorithm>
#include "problem.h"
#include "constraint_handle.h"
#include "conjunct.h"
#include "debug.h"
#include <set>
#include <boost/numeric/ublas/vector.hpp> 
#include "equalities/simpl_eq.hpp"

namespace ublas = boost::numeric::ublas; 

using namespace rfm_solver;


Problem::Problem(): 
    _normalized(true),
    _inequs(new SparseMatrix(0,0)),
    _const(new Constants(0)),
    _cnts_type(new ConstraintsType()),
    n_AllVars(0)
    {}
Problem::Problem(size_t n_AllVars1):
    _normalized(true),
    _inequs(new SparseMatrix(0,n_AllVars1)),
    _const(new Constants(0)),
    _cnts_type(new ConstraintsType()),
    n_AllVars(n_AllVars1)
    {}
Problem::Problem(const Problem & other){
	operator= (other);
}

const Problem & Problem::operator = (const Problem & other){
    if (&other==this) return *this;
    
    _normalized = other._normalized;
    _inequs.reset(new SparseMatrix(*(other._inequs)));
    _const.reset(new Constants(*(other._const)));
    _cnts_type.reset(new ConstraintsType(*(other._cnts_type)));
    n_AllVars = other.n_AllVars;
    _row_hash = other._row_hash;
    return *this;
}
 
Problem::~Problem() {}

void Problem::EliminateVariable(int var_id) {
    SparseMatrix_aptr  projective = BuildProjective(var_id);
    if (projective) {
#ifdef DEBUG_PROBLEM
        cout<<"\ninq  "<<*_inequs;
        cout<<"\nconst  "<<*_const;
        cout<<"\nproj "<<*projective;
    	cout<<"\nEliminateVariable before mult s11 "<<_inequs->size1()<<
    	" s21 "<<_inequs->size2()<<" s12 "<<projective->size1()<<
    	" s22 "<<projective->size2()<<endl<<flush;
#endif
        (*_inequs) = (*projective) * (*_inequs);
        projective->prod(*_const);  /// _inequs = projective * _inequs 
        debugf("\nproblem size(",_inequs->size1(), ", ", _inequs->size2(), ")");
        assert(_inequs->size1() == _const->size());
	_normalized = false; // who knows
    } else {
        //alredy eliminated
    }
    
    return;
};

int Problem::newEQ() {
    //check if not exceed max
    check_number_INEQs();
    //resize matrix

    _inequs->resize(nINEQs()+1, _inequs->size2()>0 ? _inequs->size2() : n_AllVars, true);
    //return number of
    _const->resize(nINEQs());
    _cnts_type->push_back(Constraint_Property(TEQ, BLACK));

    return nINEQs()-1;
};
int Problem::newGEQ() {
    //check if not exceed max
    check_number_INEQs();
    //resize matrix
    _inequs->resize(nINEQs()+1, _inequs->size2()>0 ? _inequs->size2() : n_AllVars, true);
    _const->resize(nINEQs());
    
    //return number of
    _cnts_type->push_back(Constraint_Property(TGEQ, BLACK));
    return nINEQs()-1;
};


int Problem::newGE() {
    //check if not exceed max
    check_number_INEQs();
    //resize matrix
    _inequs->resize(nINEQs()+1, _inequs->size2()>0 ? _inequs->size2() : n_AllVars, true);
    
    _const->resize(nINEQs());
    //return number of
    _cnts_type->push_back(Constraint_Property(TGE, BLACK));
    return nINEQs()-1;
};

coef_t Problem::update_coef(const size_t& row, Variable_ID v, coef_t val) {
    assert(row < _inequs->size1());
    if (v->id()+1 > n_AllVars) { // +1 since id() starts form 0 but size from 1
    	n_AllVars = v->id()+1;
    	_inequs->resize(_inequs->size1(), n_AllVars);
    }
    _normalized = false;
    return ((*_inequs)(row, v->id())+=val);
}

coef_t Problem::update_const(const size_t&row, coef_t val) {
    assert(row < _const->size());
   (*_const)[row] += val;
    return (*_const)[row];
};

coef_t Problem::get_coef(const size_t& row, Variable_ID v) const{
    assert(row < _inequs->size1());
    return (*_inequs)(row, v->id());
};
coef_t Problem::get_const(const size_t& row) const {
    assert(row < _inequs->size1());
    return (*_const)[row];
};

bool Problem::is_in_constr(const size_t& row, Variable_ID v) const{
    assert(row < _inequs->size1());
    return ((*_inequs)(row, v->id())!=0);
}


bool Problem::check_number_INEQs()const {
    if (nINEQs() > maxINEQs) {
        cerr<< "\nERROR:\n";
        cerr<< "An attempt was made to set the number of available inequality constraints to "<<nINEQs()<<".\n";
        cerr<< "The maximum number of equality constraints in a conjunction is maxINQs."<<maxINEQs<<"\n";
        cerr<< "This limit can be changed by redefining maxINQs in problem.h and recompiling.\n\n";
        exit(1);
    }
    return true;
};


//builds sparce Elimination matrix
// return true if elimination is not trivial
SparseMatrix_aptr Problem::BuildProjective(size_t column) {

#ifdef DEBUG_PROBLEM
    cout<<"\nProblem::BuildProjective: column = "<<column;
#endif
    
    assert(column<_inequs->size2());
    Coefs varCoef;
    size_t size1 = _inequs->size1();
    
    for(size_t s1=0;s1<size1; s1++)
    	varCoef.push_back((*_inequs)(s1,column));
    
    PInts ref_pos; // index of positive members of varCoef
    PInts ref_neg; // .. negative ..
    PInts ref_nul; // .. 0 .. or other than GEQ, GE

    // buld indexes
    for (size_t i = 0; i<size1;i++) {
        if (varCoef[i]>0) {
            ref_pos.push_back(i);
        } else if (varCoef[i]<0) {
            ref_neg.push_back(i);
        } else {
            ref_nul.push_back(i);
        }
    }
    //if only negative of only positive exists we
    //just delete column in matrix A
    // Also We delete this colums anyway
    // _inequs->erase_line2(column);    /// smx1
    _inequs->reset_s2(column);     /// smx2 does not have erase_line2 method
    //if we have only one type of inequality
    //we mark it as REM in _cnts_type
    // letter it will be deleted by remove_eliminated()
    
    if (ref_pos.empty() || ref_neg.empty()){
        std::for_each(
            ref_pos.begin(),
            ref_pos.end(),
            boost::bind(&Problem::reset_line,this,_1)
         );	
         std::for_each(
            ref_neg.begin(),
            ref_neg.end(),
            boost::bind(&Problem::reset_line,this,_1)
         );
        return SparseMatrix_aptr(); // empty
    }
    // reseting
    /// TODO: profile,and if too many ref_nuls
    /// move those ref_nul to other system
    /// and merge it to the main system after projection
    size_t new_size = ref_pos.size()*ref_neg.size()+ref_nul.size();
    debugf("\nprojective matrix(", new_size, ", ", size1, ") created.");
    //cout<<"\rmatrix("<< new_size<<", "<<size1<<") created."<<flush;
        
    if (new_size>=INT_PROBLEM_SIZE){
    	cerr<<"Maximum allowed internal problem size: "<<INT_PROBLEM_SIZE;
    	cerr<<"\nThe problem too complicated ("<<new_size<<", "<<size1<<")\nAborted."<<endl<<flush;
    	exit(1);
    }
    /**
    For sake of perfomance guessing of new number of
    non empty elements for the projective matrix is useful to
    implement.
   */
    SparseMatrix_aptr projective(new SparseMatrix(new_size, size1));
    ConstraintsType_aptr new_cnts_type(new ConstraintsType(new_size));	
    //buildin C for positive and negative
    // assuming that ref_pos and ref_neg at least has size 1
    int CurRow=0;
    for(size_t pp=0; pp < ref_pos.size(); pp++ ) {
        for(size_t pn = 0; pn < ref_neg.size(); pn++) {
            (*projective)(CurRow,ref_pos[pp]) = -varCoef[ref_neg[pn]];
            (*projective)(CurRow,ref_neg[pn]) = varCoef[ref_pos[pp]];
            // note. there are no equalities here
            if ((*_cnts_type)[ref_pos[pp]] == (*_cnts_type)[ref_neg[pn]]){
            	(*new_cnts_type)[CurRow] = (*_cnts_type)[ref_pos[pp]];
            } else {
            	if ((*_cnts_type)[ref_pos[pp]].type == TGE)
            		(*new_cnts_type)[CurRow] = (*_cnts_type)[ref_pos[pp]];  // > * >= substituted by >
            	else {
#ifndef NDEBUG	
            	     if ((*_cnts_type)[ref_pos[pp]].type != TGEQ){
            	     	cout<<"\nExpecting TREQ but got:"
            	     	    <<((*_cnts_type)[ref_pos[pp]].type)<<endl<<flush;
                        assert(0 && "not possible");
            	     }
#endif
            	     assert((*_cnts_type)[ref_pos[pp]].type == TGEQ);
            	     
            	    (*new_cnts_type)[CurRow] = (*_cnts_type)[ref_neg[pn]];
            	}
            }
            CurRow++;
        }
    }
    // bilding C for 0 by seting it to 1
    // assuming size of ref_nul could be 0

    for(size_t p0 = 0; p0 < ref_nul.size(); p0++) {
        (*projective)(CurRow, ref_nul[p0]) = 1;
        (*new_cnts_type)[CurRow] = (*_cnts_type)[ref_nul[p0]];
        CurRow++;
    }
    _cnts_type = new_cnts_type;
#ifdef DEBUG_PROBLEM    
    cout<<"\nbuilded c\n"<<*projective<<flush;
#endif
    return projective;
};
 
const Constraint_Property &Problem::get_ineq_type(size_t  index) const{
    assert(index <_cnts_type->size());
    return (*_cnts_type)[index];
};

Constraint_Property &Problem::ineq_type(size_t  index){
    assert(index <_cnts_type->size());
    return (*_cnts_type)[index];
};


void Problem::EliminateEquVariable(const size_t & line_id, const size_t & var_id ){
    assert(get_ineq_type(line_id).type == TEQ && "equality expected");
    
    EqSolve<SparseMatrix, Constants> eq_solver (*_inequs, *_const);
    eq_solver.Solve(line_id, var_id);
    
    reset_line(line_id);
    	
    _normalized = false;
    return;
}

size_t Problem::var_difficulty(const size_t & var_id) const {
   assert(var_id < _inequs->size2() && "invalid variable id");
    size_t pos = 0;
    size_t neg = 0;

    for(SparseMatrix::const_iterator1 it1= _inequs->begin1(var_id); it1; ++it1){
    	if (*it1<0){
    	    ++neg;
    	    continue;
    	}
    	++pos;
    }
    return pos*neg;
}

int Problem::choice_inq(const Variables &vars) const{
    // For each varible of type Wildcard_Var
    // calculate n_pos * n_neg
    // if the calulated value is less then current one
    // save the index and the value

    int index=-1;
    if (!EnableVaribleChoise)
    	return (vars.size()? 0: index);
    
    int value = std::numeric_limits<int>::max();

    for(Variables::const_iterator it= vars.begin(); it!=vars.end(); ++it){
        assert((*it)->GetType() == Wildcard_Var);
    	int c_val = var_difficulty((*it)->id());
    	if (c_val < value){
    	    index = it - vars.begin();
    	    value = c_val;
    	    if (value ==0 ) // cannot be simpler so exiting
    	    	break;  
    	}
    }
    
    return index;
}

Pair Problem::choice_eq(const Variables &vars, const PIntList_aptr & ptr2equs) const {
    
    int index=-1;
    int eq_line_id = -1;
    for(Variables::const_iterator it_var = vars.begin();it_var!= vars.end() && ptr2equs; ++it_var){ 
    	assert((*it_var)->GetType() == Wildcard_Var);
        for(PIntList::const_iterator it_eq = ptr2equs->begin();it_eq!=ptr2equs->end();++it_eq){ // EQ lines
           if (!can_eliminate(*it_eq, (*it_var)->id())) continue;
           assert(get_ineq_type(*it_eq).type == TEQ);
           index = it_var - vars.begin();
           return make_pair(index, *it_eq);
       }
    }
    return std::make_pair(index, eq_line_id);
}

bool Problem::validate_trivial(const size_t & line_id) const {
    Constraint_Property ct = (*_cnts_type)[line_id];
    const coef_t & val = (*_const)[line_id];
    switch (ct.type){
    	case TGE:{return (0 > val); }
    	break;
    	case TGEQ:{return (0 >= val); }
    	break;
    	case TEQ:{return (0 == val); }
    	break;
    	case TREM:{return true;}
    	break;
    	default:
    	assert(0 && "bad type in vector of inequalities");
    }
    return false;
};

namespace rfm_solver {

class Ptr2equBuilder{
    PIntList_aptr eq_list_;
    size_t count;
    size_t *nequs_;
public:
    Ptr2equBuilder(PIntList_aptr &eq_list, size_t * nequs):eq_list_(eq_list), count(0), nequs_(nequs) {}
    void operator ()(const Constraint_Property & ct){
          if (ct.type == TEQ){
              eq_list_->push_back(count);
              ++*nequs_;
          }
          ++count;
    }
};

}; // end namespace rfm_solver


size_t Problem::remove_eliminated(){

    debug("\nProblem::remove_eliminated before");
    prn_problem();
     
    assert(_inequs->size1() == _cnts_type->size()&& "shall by in sync before");
    assert(_const->size() == _cnts_type->size()&& "shall by in sync before");
/// physicaly erase removed lines in the matrix
    size_t nremoved = _inequs->remove_empty_s1();
    if (!nremoved) return 0;
/// remove all TREM from _cnts_type
    
    
    for(size_t i=_const->size(); i>0; --i){
    	if ((*_cnts_type)[i-1].type == TREM) {
    	    _const->erase(_const->begin()+i-1);
    	}
    }
    
    _cnts_type->erase(std::remove_if(_cnts_type->begin(), _cnts_type->end(), 
         boost::bind(std::equal_to<Constraint_Type>(),
             boost::bind(&Constraint_Property::type, _1), TREM)), _cnts_type->end());
    
    RebuildHash();
    
    debug("\nProblem::remove_eliminated after");
    prn_problem(); 
    
    assert(_const->size() == _cnts_type->size()&& "shall by in sync before");
    assert(_inequs->size1() == _cnts_type->size()&& "shall by in sync after");
    return nremoved;
}

/**
 * Now mark only lines where all coeficients is 0
 * 
 */

bool Problem::reset_trivial(){ 
    Constraint_Type ct;
    for(size_t line_id =0; line_id < _inequs->size1(); ++line_id) {
        ct = (*_cnts_type)[line_id].type;
        if(TREM == ct) continue;
        if(TINCONS == ct) return false;
   	if (is_trivial(line_id)){
   	   if(!validate_trivial(line_id)) return false;
   	   reset_line(line_id);
   	} 
   	else {
   	} 
    }
    return true;
};

inline bool Problem::is_valid(size_t line_id) const{
  Constraint_Type ct=(*_cnts_type)[line_id].type;
   return (ct!=TREM && ct!=TINCONS);
}

inline bool Problem::is_trivial(const size_t & line_id) const {
    return (_inequs->nlive(line_id)==0);
}

namespace rfm_solver {

typedef vector<Constraint_Handle> Constaraints; // Copying is cheap for Constraint_Handle

class ConstraintFill{
   Constaraints * cnts_;
public:
   ConstraintFill(Constaraints & cnts):cnts_(&cnts){}
   inline void add(Constraint_Handle & cnt){
   	cnts_->push_back(cnt);
   }
   
};


class RedundantMarker{
    Problem_ptr problem_;
    Constraint_Handle cnt_;
    size_t removed_count_;
    bool was_incons;
    bool red_only_;
    bool was_new_eq;
    Conjunct_aptr in_gist_; /// inconsistent RED conjunct putted here
    Conjunct_ptr parent_conj_; /// inconsistent RED conjunct putted here
    
public:
    
    RedundantMarker(Problem_ptr p, Constraint_Handle & c, bool ro, Conjunct_ptr parent):
        problem_(p), cnt_(c), removed_count_(0), was_incons(false), red_only_(ro), was_new_eq(false), parent_conj_(parent){}
    
    /// Save black from being redundant
    int decide_color(const Constraint_Handle & c0, const Constraint_Handle & c1, int redundant_line);
    
    bool mark (const Constraint_Handle & c, int (Problem::*choise_red)(const size_t &, const size_t &) const, bool (Constraint_Handle::*coef_eq)(const Constraint_Handle &)const );
    
    bool process_redundant(const Constraint_Handle  * c2, int &);
    
    size_t n_removed()const {return removed_count_;}
    void reset_constraint(Constraint_Handle & c){
    	cnt_ =  c;
    }
    bool consistent()const {return !was_incons;}
    bool new_eq()const {return was_new_eq;}
    Conjunct_aptr in_gist(){return in_gist_;}
    void prn_debug(const Constraint_Handle * c1, const Constraint_Handle * c2);
    		
};
}; /// end of namespace

/**
decide_color
tunes made desision according to constrasint colors
rline 	0	zero deleted	
	BLACK	RED	0
BLACK	0	0	
RED	-1	0	
1			
rline 	1	one deleted	
	BLACK	RED	0
BLACK	1	-1	
RED	1	1	
1			
rline 	2	don't matter	
	BLACK	RED	0
BLACK	2	0	
RED	1	2	
rline 	2	don't matter	
*/
int RedundantMarker::decide_color(const Constraint_Handle & c0, const Constraint_Handle & c1, int redundant_line) {

    Constraint_Color color0  = c0.get_type().color;
    Constraint_Color color1  = c1.get_type().color;
    if (redundant_line == FIRST) {
        if (color0 == BLACK && color1 == RED)
            return NOT_REDUNDANT;  // black cannot be redundant
        return FIRST;
    } else if (redundant_line == SECOND) {
        if 	(color0 == RED && color1 == BLACK)
            return NOT_REDUNDANT; // black cannot be redundant
        return SECOND;
    } else if (redundant_line == DONT_CARE ) { // don't care
        if (color0 == RED && color1 == BLACK)
            return FIRST; // choose red
        if (color0 == BLACK && color1 == RED)
            return SECOND;  // choose red
    } else if (redundant_line == INCONS) { 
    /// thouse inconsistent who definitely will be in gist removed from other RED constraints
        if (color0==RED || color1==RED) {
    	    if(!in_gist_) in_gist_ = parent_conj_->create();
    	    if (color0==RED && color1==BLACK) {
    	    	in_gist_->add_CNT(c0); /// TODO:: may be implement swap of inequalities
    	    	return FIRST;
    	    } else if (color0==BLACK && color1==RED) {
    	    	in_gist_->add_CNT(c1);
    	    	return SECOND;
    	    }
    	    assert(0 && "2 READ int gist are inconsistent");
    	    return INCONS;
    	} else {
    	    assert(0 && "BLACK int gist are inconsistent");
    	}
        return INCONS;
    }
    return redundant_line;
}

void RedundantMarker::prn_debug(const Constraint_Handle * c1, const Constraint_Handle * c2) {
#ifdef DEBUG_PROBLEM
    PrettyPrintVisitor prn;
    cout<<"\n";
    if (c1) {
        Constraint_Handle  cc1 = *c1;  /// because of const
        cc1.Accept(prn);
        cout<<"("<<cc1.get_index()<<")";
    } else {
        cout<<" NULL";
    }
    cout<<" == ";
    if (c2) {
        Constraint_Handle  cc2 = *c2;
        cc2.Accept(prn);
        cout<<"("<<cc2.get_index()<<")";
    } else {
        cout<<" NULL";
    }
    cout<<endl<<flush;
#endif
}

bool RedundantMarker::mark(const Constraint_Handle & c2, 
    int (Problem::*choise_red)(const size_t &, const size_t &) const, 
    bool (Constraint_Handle::*coef_eq)(const Constraint_Handle &)const  ){
#ifdef DEBUG_PROBLEM
    cout<<"\nmark input:";
    prn_debug(&cnt_, &c2);
#endif
    
if (cnt_.empty()) {
        cnt_ = c2;
        return true;
    }
    assert(cnt_.get_index() != c2.get_index() && "trivial redundancy");
    
    int redundant_line = NOT_REDUNDANT;
    if ((cnt_.*coef_eq)(c2)) {
        if (red_only_) {
            if (NOT_REDUNDANT!=(redundant_line=(problem_->*choise_red)(cnt_.get_index(),c2.get_index())))
                redundant_line = decide_color(cnt_, c2, redundant_line);
        } else {
            redundant_line=(problem_->*choise_red)(cnt_.get_index(),c2.get_index());
        }
    }
    if (NOT_REDUNDANT != redundant_line) {
        if (!process_redundant(&c2, redundant_line)){ // inconsistency found
            was_incons = true;
            return false;

        }

#ifdef DEBUG_PROBLEM
    cout<<"\nafter process_redundant:";
    prn_debug(&cnt_, &c2);
#endif

    } else { // not redundant, just move next
#ifdef DEBUG_PROBLEM
        cout<<"\n=== different ===";
#endif
        cnt_ = c2;
    }
    return true;
}

bool RedundantMarker::process_redundant(const Constraint_Handle  * c2, int& redundant_line){
#ifdef DEBUG_PROBLEM
    cout<<"\nprocess_redundant examined lines: ("<<cnt_.get_index()<<","<< c2->get_index()<<")";
    cout<<"\n===redundant or incons line: "<<redundant_line<<endl<<flush;
#endif
    bool consistent = true;
    switch (redundant_line) {
    case FIRST: {
            problem_->reset_line(cnt_.get_index());
            ++removed_count_;
            cnt_ = *c2;
            break;
        }
    case SECOND:
    case DONT_CARE: {
            problem_->reset_line(c2->get_index());
            ++removed_count_;
            
            break;
        }
    case INCONS: {
            problem_->change_line_type(cnt_.get_index(), TINCONS);
            cnt_.clear();
            problem_->change_line_type(c2->get_index(), TINCONS);
            consistent = false;
            
            break;
        }
    case DK2EQU: {  // current to equality, reset latest
            problem_->change_line_type(cnt_.get_index(), TEQ);
            problem_->reset_line(c2->get_index());
            ++removed_count_;
            was_new_eq = true;
            
            break;
        }
    case BOTH: {   /// for gist only
            problem_->change_line_type(cnt_.get_index(), TREM);
            cnt_.clear();
            problem_->change_line_type(c2->get_index(), TREM);
            
            break;
    }
    default:
        assert(0 && "unexpected value received from choice_redundant");

    }
#ifdef DEBUG_PROBLEM
    if (consistent)
    	cout<<"\nline consistent"<<endl;
    else
    	cout<<"\nline inconsistent"<<endl;
#endif

    return consistent;
}

RedundantPOut Problem::eliminate_p_redundant(Conjunct_ptr c, bool red_only){
    RedundantPOut status = 
        eliminate_p_same(c, red_only);
    
    if (status.consistent){
    	RedundantPOut ops_status = eliminate_p_opposite(c, red_only);
    	if (!ops_status.consistent) {
    	    status.update(ops_status);
    	    return status;
    	} 
    	status.update(ops_status);
    }
    
    if (status.nremoved > 0)	/// some redundancy was found
    	remove_eliminated();
    return status;
};

/** 
  * x > a x = b, (a < b)? ( x=b : false )
   * x >=a x = b, (a <=b)?( x=b : false)
  * x = a x = b, false 
  **/
  
int Problem::choice_equ_inequ(const int & eq_no, const coef_t &eq_val, const coef_t &inq_val, const Constraint_Type & inq_type)const {
    if (inq_type == TEQ && inq_val != eq_val)
   		return INCONS;
    if (inq_type == TGEQ) 
       return ((inq_val<= eq_val) ? eq_no:INCONS);
    if (inq_type == TGE) 
	return ((inq_val < eq_val) ? eq_no:INCONS);
    assert(0 && "unexpected type");
    return INCONS;
}

/** choise a nonredundant line
 * x > val0 & x > val1 if (val0 less than val1) --> #1 
 * x==val0 & y==val1 -->unknown
 * x==val x>=val -->  
 */
int Problem::choice_redundant(const size_t & line_id0, const size_t & line_id1 )const{
   const coef_t & val0 = get_const(line_id0);
   const coef_t & val1 = get_const(line_id1);
   Constraint_Property ct0 = (*_cnts_type)[line_id0];
   Constraint_Property ct1 = (*_cnts_type)[line_id1];
   if (ct0.type==TREM || ct1.type==TREM) 
   	return -1;
   if (val0!=val1){
   	if (ct0.type==TEQ )
   	    return (choice_equ_inequ(1, val0, val1, ct1.type));
   	if (ct1.type==TEQ )
   	    return (choice_equ_inequ(0, val1, val0, ct0.type));	
   	return (val0>val1)? 1:0 ; 
   }
   // calculate using the table
   int res = redundant_tbl[ct1.type][ct0.type];
    return res;
}
/** choise a nonredundant line
  *  b == -a 
  * -x >=a & x >= b  --> b <= x <= -a --> (x = b)
  * -x >=a & x = b  --> x=b & x<=-a --> (x = b)
  * -x = a & x = b	 --> (x = b)
  *         other cases --> false
  * Assuming a != b:
  * -x >(=) a & x >(=) b  --> b <= x (=)< -a --> if -a > b --> false else not_redundant
  * x >(=) a & -x >(=) b  --> a <= x (=)< -b --> if (-b > a) == ( -a > b ) --> false else not_redundant
  * -x = a & x >(=) b or
  *  x = a & -x >(=) b    --> if -a > b leave the equality else false
  *  x = a & -x = b  -->  if (a == -b) leave the equality else false
  * TODO: figure out all of the theory
  */

int Problem::choice_redundant_opposite(const size_t & line_id0, const size_t & line_id1 ) const {
   
   const coef_t & val0 = get_const(line_id0); // the only difference from choice_redundant
   const coef_t & val1 = get_const(line_id1);
   Constraint_Property ct0 = (*_cnts_type)[line_id0]; 
   Constraint_Property ct1 = (*_cnts_type)[line_id1];
   if (ct0.type==TREM || ct1.type==TREM) 
   	return NOT_REDUNDANT;
   if (val0!=-val1){
   	if (ct0.type==TEQ )
   	    // equ 
   	    return (choice_equ_inequ(1, -1*val0 /*eq*/, val1, ct1.type)); //implements -eq const > inq const
   	if (ct1.type==TEQ )
   	    return (choice_equ_inequ(0, -1*val1, val0, ct0.type));	
   	return (-val0>val1)? NOT_REDUNDANT:INCONS; 
   }
   // calculate using the table
   int res = opposite_redundant_tbl[ct1.type][ct0.type];
    return res;
}

int Problem::normalize_by_var(size_t var_id, Marker & nm){
    // deviding every variable for the coeficient of the variable
    // to be eliminated
    size_t line_id;
    coef_t local_k;
    for(SparseMatrix::iterator1 it1=_inequs->begin1(var_id); it1; ++it1){
    	line_id = it1.get_s1();
    	if (nm.IsMarked(line_id)) continue; // inequality already normlized
    	nm.Mark(line_id, var_id);
    	if (get_ineq_type(line_id).type== TEQ) {
    	    local_k = *it1; // to eliminate case -x = a
    	} else {
        	local_k = (*it1<0)? -1*(*it1):*it1; // got non zero coefficient, take modulo
        }
        for(SparseMatrix::iterator2 it2=_inequs->begin2(it1.get_s1()); it2; ++it2){ 
    	    *it2 = *it2/local_k;
        }
        (*_const)[line_id] /= local_k;
    }
    
    return 0;
}

/**
    int Problem::normalize_problem(Variables &);
    create int vector with the same size as number of inequalities
    fill it with 0 which means all eventualities are not normalized 
    for each columns (live variable) normalize all rows (inequalities)
    if it is not marked by 1 in the vector.
    mark it as normalized by assign corresponding element of the vector to 
    Variable.id() value. The value is always more than 0 and equal to the column number.
    Always return true;
**/

/// TODO: refactor it. Use row ID instead of variables
int Problem::normalize_problem(const Variables & locals){
    Marker nm(size());
    size_t index = 0;
    while(!nm.IsAllMarked() && index < locals.size()){
    	normalize_by_var(locals[index]->id(), nm);
    	++index;
    }
    RebuildHash();
    _normalized = true;
    return 0;
}

void Problem::RebuildHash(){
    _row_hash.clear();
    size_t nrows = this->size();
    for(size_t irow=0; irow < nrows;++irow){
    	_row_hash.insert(RowHash::value_type(irow, _inequs->coefs_hash(irow)));
    }
}

bool Problem::HashRehashRow(size_t  line_id){
    Row_index& row_index = _row_hash.get<first>();
    Row_index::iterator it = row_index.find(line_id);
    if (it!= row_index.end()) {
    	Row_index::value_type curr = *it;
    	curr.second= _inequs->coefs_hash(line_id);
        row_index.replace(it, curr);
        return true;
    }
    return false;
}

bool Problem::HashRemoveRow(const size_t& line_id){
    Row_index::iterator it=_row_hash.get<first>().find(line_id);
    if (it != _row_hash.get<first>().end()){
        _row_hash.erase(it);
        return true;
    }
    return false;
}

std::pair<bool, size_t> Problem::find_opposite_row(size_t line_id) const {
    size_t opposite_hash = _inequs->coefs_hash_opposite(line_id);
    Hash_index::const_iterator it=_row_hash.get<second>().find(opposite_hash);
    if (it != _row_hash.get<second>().end()){
    	return std::make_pair(true, it->first);
    }
    return std::make_pair(false, 0);
}

void Problem::reset_line(const size_t & line_id){ 

    if (get_ineq_type(line_id).type != TREM){
        change_line_type(line_id, TREM);
    }
    _inequs->reset_s1(line_id);
    /// cannot remove values from hash so lets rehash it
    HashRehashRow(line_id);
}

/**
******** Marker implementation *******************
**/
Marker::Marker(size_t s):_count(0), _vnorm(s,0){};

bool Marker::IsAllMarked() const {
	return (_vnorm.size() == _count)? true:false;
}

void Marker::Reset(size_t s){
    _count = 0;
    _vnorm.resize(s);
    _vnorm.assign(_vnorm.size(), 0);
};

void Marker::Mark(size_t index, size_t val){
    assert(index < _vnorm.size());
    if (val == 0 && _vnorm[index] != 0) --_count;	
    else if (val != 0 && _vnorm[index] == 0) ++_count;
    else {/* nothing */}
    _vnorm[index] = val;
    return;
}

inline const size_t & Marker::operator()(size_t index) const{
    assert(index< _vnorm.size());	
    return _vnorm[index];
};

inline bool Marker::IsMarked(size_t index) const {
     assert(index < _vnorm.size());
     return (_vnorm[index] == 0) ? false:true;
}

void Problem::change_line_type(size_t line_id, Constraint_Type new_type) {
    Constraint_Type & old_type = ineq_type(line_id).type;
    if (old_type == new_type){
        return;
    }
     old_type = new_type;
}

bool Problem::operator == (const Problem & other) const {
     if (this != &other){
     	return (
     	  _normalized == other._normalized &&
     	  *_cnts_type == *other._cnts_type &&
     	  *_inequs == *other._inequs
     	);
     }
     return true;
 }

bool Problem::operator < (const Problem & other) const{
    if (this != &other){
     	return (
     	nINEQs() < other.nINEQs() ||
 
     	*_inequs < *other._inequs ||
     	*_cnts_type < *other._cnts_type 
   
     	);
     }
     return false;
}

void rfm_solver::Problem::erase_line( size_t line_id )
{
    _inequs->erase_s1(line_id);
    _cnts_type->erase(_cnts_type->begin()+line_id);
    assert(_inequs->size1()==_cnts_type->size());
}

Problem::Ptr2Eq rfm_solver::Problem::create_ptr2equs( ) const
{
    PIntList_aptr ptr2eq(new PIntList);
    size_t nequs = 0;
    Ptr2equBuilder ptr2eq_builder(ptr2eq, &nequs);
    for_all(*_cnts_type, ptr2eq_builder);
    return make_pair(ptr2eq, nequs);
}

RedundantPOut Problem::eliminate_p_same(Conjunct_ptr c, bool red_only) {
    
    /** 
    If size of hash less or equal 1 retur true
    Get element from the hash, cunstruct Constraint_Handle
    pass to RedundantMarker 
    for all next elements, apply mark method of the RedundantMarker
    pass to redunadant marker only potentialy equal values
    **/
    
    if ( this->size() <= 1 )
        return RedundantPOut(true,0);
    
    Constraint_Handle c1;
    Constraint_Handle c2;

    if (!_normalized)
        normalize_problem(c->all_variables());
    
    assert(this->size() == _row_hash.size()&& "mast be in sync");

    
    /// create  a view indexed by hash
    
    const Hash_index& hash_id=_row_hash.get<second>();
    bool consistent = true;
     

    Hash_index::const_iterator it = hash_id.begin();
    Hash_index::const_iterator hash_end = hash_id.end();
    size_t h2, h1 = it->second; /// second is hash
    
    c1 = Constraint_Handle(c, it->first); /// first is row #
    RedundantMarker r_mrk(this, c1, red_only, c);
    ++it;

    while(it!=hash_end && consistent){
       h2 = it->second;
       if (h2!=h1) {
           c1 = Constraint_Handle(c, it->first);
           r_mrk.reset_constraint(c1);
           ++it;
       } else {
       c2 = Constraint_Handle(c, it->first);
        ++it; // get next iterator because current may be invalidated by the next line
       
       consistent = r_mrk.mark(c2, &Problem::choice_redundant,
            &Constraint_Handle::coeficients_equal);
      }
       h1 = h2;
    }


#ifdef DEBUG_PROBLEM
    cout<<"\neliminate_p_same return:"<<(r_mrk.consistent()? "TRUE":"FALSE")<<endl<<flush;
    prn_problem();
    if (r_mrk.new_eq()){
    	cout<<"\neliminate_p_same: new equality created";
    }
    if (r_mrk.in_gist()) {
    	cout<<"\neliminate_p_same: constraints in gist:";
    	debugg(r_mrk.in_gist());
    }
#endif

    return RedundantPOut(r_mrk.consistent(), r_mrk.n_removed(), r_mrk.new_eq(), r_mrk.in_gist());
}

RedundantPOut Problem::eliminate_p_opposite(Conjunct_ptr c, bool red_only) {

    if ( this->size() <= 1 )
        return RedundantPOut(true,0);

    // normalize by myVars

    if (!_normalized)
        normalize_problem(c->all_variables());   // bugfix
    assert(this->size() == _row_hash.size()&& "mast be in sync");
    

    typedef std::set<size_t, std::less<size_t> > RowsSet;
    RowsSet done; /// set of alredy visited lines (rows)
    
    Constraint_Handle c1, c2;
    RedundantMarker r_mrk(this, c1, red_only, c);
    bool consistent = true;
    
    const Row_index& row_id=_row_hash.get<first>();
    Row_index::const_iterator it = row_id.begin();
    Row_index::const_iterator row_end = row_id.end();
    
    while(it!=row_end && consistent){
    	/// find opposite row for current 
    	/// if it's not in the done set
    	///    call simplification for it
    	///   put both to the done set
    	size_t r1 = it->first;
    	++it;
    	if (done.find(r1)!=done.end() && is_valid(r1))
    	    continue;
    	pair<bool, size_t> res = find_opposite_row(r1);
    	size_t r2 = res.second;
    	if ( res.first && is_valid(r2)){
    	   
    	    c1=Constraint_Handle(c, r1);
    	    c2=Constraint_Handle(c, r2);
    	    r_mrk.reset_constraint(c1);
    	    consistent = r_mrk.mark(c2, &Problem::choice_redundant_opposite,
    	        &Constraint_Handle::coeficients_opposite_equal);
    	    done.insert(r2);
    	} else if (res.first){ /// found opposite is alredy removed
    	    done.insert(r2);
    	};
    	done.insert(r1);
    }

#ifdef DEBUG_PROBLEM
    cout<<"\neliminate_p_opposite return:"<<(r_mrk.consistent()? "TRUE":"FALSE")<<endl<<flush;
    prn_problem();
    if (r_mrk.new_eq()){
    	cout<<"\neliminate_p_opposite: new equality created";
    }
    if (r_mrk.in_gist()) {
    	cout<<"\neliminate_p_opposite: constraints in gist:";
    	debugg(r_mrk.in_gist());
    }
#endif

    return RedundantPOut(r_mrk.consistent(), r_mrk.n_removed(), r_mrk.new_eq(), r_mrk.in_gist() );
};

void RedundantPOut::update(const RedundantPOut & o) {
    consistent = consistent && o.consistent;
    nremoved += o.nremoved;
    eq_created = eq_created || o.eq_created;
    if (in_gist && o.in_gist) {
        in_gist = in_gist->merge_conjs(in_gist, o.in_gist);
    } else if (o.in_gist) {
        in_gist = o.in_gist;
    }
}

bool Problem::solve_eqs() {
	using namespace boost::numeric::ublas;
    ublas::matrix<coef_t> A;
    ublas::vector<coef_t> c;
	EqMapingData md;
	if (export_eq(A, c, md)){
	//	std::cout<<"\nExtracted A = "<<A;
	//	std::cout<<"\nExtracted c = "<<c<<std::endl;
		if(solve_equ(A, c)){
			return import_eq(c, md);
		}
	}
	return false;
}

    
    template <class MX, class V> 
	bool Problem::export_eq(MX & mx, V & c, EqMapingData & md) const {
		size_t n_equs = nEQUs();
		size_t n_vars = _inequs->size2();
		if (n_equs == 0)
			return false;

		c.resize(n_equs);
		EqMapingData::VarMap_aptr vm(new EqMapingData::VarMap);
		EqMapingData::EqMap_aptr eqm(new EqMapingData::EqMap);
		std::map<size_t, size_t>  var2col;
		typedef std::map<size_t, size_t>::iterator It; 
		size_t eq_count = 0;
		size_t var_count = 0;
		size_t all_cnts = _cnts_type->size();
		for(size_t id = 0; id < all_cnts; ++id){
			if((*_cnts_type)[id].type != TEQ)
				continue;
			eqm->push_back(id);
			for(size_t prob_var = 0; prob_var < n_vars; ++prob_var){
				coef_t val = (*_inequs)(id, prob_var);
				if (val != 0){
					It it = var2col.find(prob_var);
					if (it==var2col.end()){
						var2col[prob_var] = var_count;
						mx.resize(n_equs, var_count+1);
						vm->push_back(var_count);
						assert(vm->size() == var_count+1);
						mx(eq_count, var_count) = val;
						++var_count;
					} else {
						mx(eq_count, it->second) = val;
					}
				}
			}
			if (var_count > n_equs)
				return false;		//canot solve such system
			c(eq_count) = (*_const)[id];
			++eq_count;
		}
		md.var_map(vm);
		md.eq_map(eqm);	
		assert(eq_count == n_equs);
		return true;
	}
	template <class V>
	bool Problem::import_eq(const V & c, const EqMapingData & md) {
		if (c.empty())
			return false;
		EqMapingData::VarMap_aptr vm = md.var_map() ;
		EqMapingData::EqMap_aptr eqm = md.eq_map();
		/// forall equlities id
		size_t count = 0;
		assert(c.size() == vm->size());
		for(EqMapingData::EqMap::const_iterator it_eq=eqm->begin(); it_eq!=eqm->end();++it_eq, ++count ){
			/// write zero to all line
			for(size_t var = 0; var<_inequs->size2(); ++var)
				(*_inequs)(*it_eq, var) = 0;
			/// get index of the next variable
			/// put 1 there
			(*_inequs)(*it_eq, (*vm)[count]) = 1;
			(*_const)[*it_eq] = c(count);
			HashRehashRow(*it_eq);
		}
		assert(count = c.size());
		return true;
	}
